import sys
import os
import pickle
from processing import Process
import time
import random
import string

def para_process (arg,para_dict,para_type,message):
#This function is used to get input parameters from user
	if len(arg) == 0:
		print >> sys.stderr, "Input error!" + message
		return -1
	try:
		for i in range(0,len(arg),2):
			k1 = arg[i]
			k2 = arg[i+1]
			if para_type[k1] in ['i','int']:
				para_dict[k1] = int (k2)
			elif para_type[k1] in ['f','float']:
				para_dict[k1] = float (k2)
			else: 
				para_dict[k1] = k2

		for k in para_dict.keys():
			if para_dict[k] == '' and para_type[k] != '':
				print >> sys.stderr, "Input error!\nPlease input " + k 
				return -1				
	except:
		print >> sys.stderr, "Input error!" + message
		return -1
	return 1 

def printError(message): 
	print >> sys.stderr, message
	return -1

def chromList(M=1):
	chrlist = []
	for i in range (1, 25):
		chr = "chr" + str(i)
		if i == 23: chr = "chrX"
		elif i == 24: chr = "chrY"		
		chrlist.append(chr)
	##updated on 04-Dec-10
	if M==1:
		chrlist.append("chrM")
		
	return chrlist

## Wrote on 28-May-2012 to perform spearman correlation test in R
def spearman_test(x,y):
	N = 50
	common_file = ''.join(random.choice(string.ascii_uppercase + string.digits) for x in range(N))
	x_file = common_file + "__X"
	x_out = open(x_file,'w')

	y_file = common_file + "__Y" 
	y_out = open(y_file,'w')
	
	z_file = common_file + "__Z"  #''.join(random.choice(string.ascii_uppercase + string.digits) for x in range(N))
	
	for i in range(len(x)):
		x_out.write(str(x[i]) + "\n")
		y_out.write(str(y[i]) + "\n")
	x_out.close()
	y_out.close()
	
	cmd = "spearman.R %s %s %s" %(x_file, y_file, z_file)
	os.system(cmd)
	
	ct = open(z_file).readline().split()
	
	#rmcmd = "rm %s %s %s" %(x_file, y_file, z_file)
	#os.system(rmcmd)
	
	os.remove(x_file)
	os.remove(y_file)
	os.remove(z_file)

	return [ float(ct[0]), float(ct[1]) ]

## Spearman partial correlation test
def spearman_ptest(x,y,z):
	N = 50
	common_file = ''.join(random.choice(string.ascii_uppercase + string.digits) for x in range(N))
	x_file = common_file + "__X"
	x_out = open(x_file,'w')

	y_file = common_file + "__Y"  #''.join(random.choice(string.ascii_uppercase + string.digits) for x in range(N))
	y_out = open(y_file,'w')
	
	z_file = common_file + "__Z"  #''.join(random.choice(string.ascii_uppercase + string.digits) for x in range(N))
	z_out  = open(z_file,'w')

	o_file = common_file + "__O" #''.join(random.choice(string.ascii_uppercase + string.digits) for x in range(N))

	for i in range(len(x)):
		x_out.write(str(x[i]) + "\n")
		y_out.write(str(y[i]) + "\n")
		z_out.write(str(z[i]) + "\n")
	x_out.close()
	y_out.close()
	z_out.close()
	
	cmd = "partial_spearman.R %s %s %s %s" %(x_file, y_file, z_file, o_file)
	os.system(cmd)
	
	ct = open(o_file).readline().split()
	
	#rmcmd = "rm %s %s %s %s" %(x_file, y_file, z_file, o_file)
	#os.system(rmcmd)
	
	os.remove(x_file)
	os.remove(y_file)
	os.remove(z_file)
	os.remove(o_file)

	return [ float(ct[0]), float(ct[1]) ]


def chromList_new(l = range(1,26), type = "UCSC"):
## Type: UCSC, Ensembl
	chrlist = []
	if type.lower()[0] == 'u': text = "chr"
	elif type.lower()[0] == 'e': text = ""
	else:
		print >> sys.stderr, "Type not recognized in chromList_new()"
		sys.exit()
		
	for i in l:
		if i == 25:
			if text == "chr": chr = "chrM"
			elif text == "": chr = "MT"		
		else:
			n = str(i)
			if i == 23: n = "X"
			elif i == 24: n = "Y"
			chr = text + n
		chrlist.append(chr)
	return chrlist


def get_chrominfo(path="/data2/thong/human_index/chromInfo_hg19.txt", subdict = 0):
	dict = {}
	lines = getStream(path)
	chrlist = chromList_new()
	for l in lines:
		items = l.split("\t")
		chr = items[0]
		length = int(items[1])
		if chrlist.count(chr) == 0: continue
		if subdict == 0:
			dict[chr] = length
		else:
			dict[subdict[chr]] = length
	return dict

## This func get number of reads overlapping each base in the human genome
def get_base_dict(bedfile, sep="\t", skip= 0):
	dict = {}
	lines = getStream(bedfile, skip)
		
	for line in lines:
		items = line.split(sep)
		chr   = items[0]
		if not dict.has_key(chr):
			dict[chr] = {}
		start = int(items[1])
		end   = int(items[2])
		for i in range( start, (end + 1) ):
			if not dict[chr].has_key(i): dict[chr][i] = 1
			else: dict[chr][i] += 1
		
	return dict

def getStream(path, header=0):
	stream = open(path)
	for i in range(header):
		stream.readline() #skip header
	return stream.readlines()
	stream.close()

def GC (seq):
	a = seq.upper().count("A")
	t = seq.upper().count("T")
	g = seq.upper().count("G")
	c = seq.upper().count("C")
	try:
		return 100 * ( (g+c)/(a+t+g+c+0.0) )
	except:
		return ""
		

def ambiguity (seq):
	a = seq.upper().count("A")
	t = seq.upper().count("T")
	g = seq.upper().count("G")
	c = seq.upper().count("C")
	l = len(seq.strip())
	return 100 * ( (l - a -t - g - c)/(l+0.0) )	
	
#def nucleBlock (fdict, rdict, start, end, size):
## This divides a region from start -> end to many equal block (size),
## then calculates the nucleosome occupancy in each block 
#	flist = []
#	rlist = []
#	numblock =int( round( (end - start)/float(size) ) )
#	for i in range (numblock):
#		x1 = start + i * size
#		x2 = start + (i+1) * size
#		[fnum, rnum] = countReadsBoth(x1, x2, fdict,rdict)
#		flist.append( fnum )
#		rlist.append( rnum )
#	return [flist, rlist]


def nucleBlock (fdict, rdict, start, end, size, slide = -1):
## Use a sliding window with length = size and sliding step = slide 
## to scan the region from start to end the occupancy of nucleosome

	if slide == -1: slide = size 

	flist = []
	rlist = []
	for i in range (start, end, slide):
		x1 = i			# start position of sliding window
		x2 = x1 + size	# end   position of sliding window
		[fnum, rnum] = countReadsBoth(x1, x2, fdict,rdict)
		flist.append( fnum )
		rlist.append( rnum )
	return [flist, rlist]


## Count #reads belonging to each equal-size window cross genome.
## shift = 1/2 nucleosome size. skip: number of top lines to skip
def countWindow(chromfile, bedfile, winsize, shift = 75, skip = 0, sep = '\t'):
	enrich = {}

	## Get chromosome sizes based on UCSC chromsome size information
	chrinfo = open(chromfile).readlines()
	chrdict = {}
	chroms = chromList()
	for l in chrinfo:
		arr = l.split(sep)
		chr = arr[0]
		if chroms.count(chr) == 0: continue #not in the chrom list
		size = int(arr[1])
		chrdict[chr] = size


	## Get bed dict
	print >> sys.stderr, "Processing each tag from BED...."
	bed = getStream(bedfile, header = skip)

	#process each read (tag) in bed file
	for tag in bed:
		temp = tag[0:-1].split(sep) #remove \n and split
		chr = temp[0] #chromsome name in the first column of BED file
		
		if chroms.count(chr) == 0: continue #not in the chrom list 1, 2,..., 22, X, Y, M
		
		#initiation chr list (by a zeros vector) that stores the number of reads in each window of chr
		#the number of reads in a window is stored as an item of chr list
		if not enrich.has_key(chr):			
			numwin = chrdict[chr]/winsize + 1  #the last window may not have the size of 'winsize'					
			enrich[chr] = [ 0 for i in range(numwin) ]

		#find the center of nucleosome by moving 75bp forward or backward
		if temp[5] == '+': #forward strand
			center = int(temp[1]) + shift
		elif temp[5] == '-': #reverse strand
			center = int(temp[2]) - shift
		else: 
			print >> sys.stderr, "Something wrong with strand!"
			sys.exit()
			
		win = center / winsize #get the index (window) of list that the tag belongs to
		
		## added on 5-Nov to remove any reads not belonging to the list of windows of current chromsome
		if win < len(enrich[chr]):			
			enrich[chr][win] += 1

	return enrich


def countReads(x1, x2, dict):
#This is to count number of reads mapped to a specific region x1 -> x2
	val = 0
	for i in range (x1, x2):
		if dict.has_key(i):
			val += dict[i]
	return val

def countReadsBoth(x1, x2, fdict,rdict):
#This is to count number of reads (forward or reverse) mapped to a specific region x1 -> x2
	fval = 0
	rval = 0
	for i in range (x1, x2):
		if fdict.has_key(i):
			fval += fdict[i]
		if rdict.has_key(i):
			rval += rdict[i]
	return [fval, rval]

def wigDict(dict, path, chr, header = 1):
#This is to get a dictionary of nucleosome score for a chromosome
#key: chromosome; value: list of lines of the chromosome
	dict[chr] = getStream(path,header)

def changeName(name):
	name = name.replace ("(",'')
	name = name.replace (")",'')
	name = name.replace ("/",'_')
	name = name.replace (" ",'_')
	return name

def bedDict (dict,path, col=0, change_name = False, header = 1, sep = '\t'):
#This is to get a chromosome dictionary for a bed file
#col: is the column used for key value of dict
	lines = getStream(path,header)
	if len(lines) == 0: 
		return False
	for line in lines:
		key = line[0:-1].split(sep)[col]
		##key = key.replace('\n','')
		if change_name: 
			key = changeName(key)
		if not dict.has_key(key):
			dict[key] = []
		dict[key].append(line)

def writeList( arr, file, endline = '\n', mode = 'w' ):
#Write a list into file
	out = open(file,mode)
	for item in arr:
		out.write( str(item) + endline )	
	out.close()		

def listDir(path,pattern = '', ext=''):
#List directory content; mode is file or dir; pattern of the filename, ext is extention

	arr = os.listdir(path)
	results = []
	for item in arr:
		if ext != '':
			if item[-(len(ext)):] == ext:
				results.append(item)
		elif item.find (pattern)!= -1:
			results.append(item)
	results.sort()
	return results

def strandDict(nudict, col, scol):
#col is the column of start
#scol is the column of strand
	fdict = {} #forward dict
	rdict = {} #reverse dict
	for chr in nudict.keys():
		fdict[chr] = {}
		rdict[chr] = {}
		for line in nudict[chr]:
			temp = line.split()
			start = int(temp[col])
			end = int(temp[col+1])
			strand = temp[scol]
			
			if strand == "+":
				if not fdict[chr].has_key(start): fdict[chr][start] = 1
				else: fdict[chr][start] += 1
			elif strand == "-":
				if not rdict[chr].has_key(end): rdict[chr][end] = 1
				else: rdict[chr][end] += 1
			else:
				print >>sys.stderr, "strandDict func: BED file error !"
				sys.exit()
	
	return [fdict,rdict]




################################################
## Get the dictionary from BED file in this form
## {'chr1':{x1:y1, x2:y2,...},... }
## xi is the center position of nucleosome
## yi is number of reads in xi
## NTT on 08-Sep-2010
################################################
def twostrandDict(bed, col=1, scol=5, shift = 75, sep = "\t", read_weight=False):
#col is the column of start
#scol is the column of strand
#shift from start or end to center = 1/2 of nucleosome length = 75
#read_weight: added on 14-Jan-2012 to incorporate multi-reads
	nudict =  nuDict(bed)
	outdict = {}
	
	for chr in nudict.keys():
		outdict[chr] = {}
		for line in nudict[chr]:
			temp = line[0:-1].split(sep)
			start = int(temp[col])
			end = int(temp[col+1])
			strand = temp[scol]
			
			##===added on 14-Jan-2012 to include multi-reads
			weight = 1
			if (read_weight == True):
				weight = float(temp[4]) ## the weight of a read ranging from 0 to 1
			##===END added
			
			if shift == -1: shift = abs( int( (end - start)/2 ) )
			if strand == "+":
				center = start + shift
			elif strand == "-":
				center = end - shift
			else:
				print >>sys.stderr, "strandDict func: BED strand error !"
				sys.exit()
			
			## Added on 6-Sep-2011 ######
			if center <=0: center = 1
			## END Added ################
			
			if not outdict[chr].has_key(center): outdict[chr][center] = weight
			else: outdict[chr][center] += weight

	return outdict



################################################
## Get the dictionary from BED file in this form
## {'chr1':{x1:y1, x2:y2,...},... }
## xi is the position in the genome
## yi is number of tags overlapping xi
## NTT on 10-August-2012
################################################
def tag_overlap(bed, extend = 200, read_weight=False):
#extend from 5'end of read to create tag
#read_weight: added on 14-Jan-2012 to incorporate multi-reads
	outdict = {}
	totalcount = 0
	for l in open(bed):
		totalcount += 1
		it = l.strip().split("\t")
		chr = it[0]
		if not outdict.has_key(chr): outdict[chr] = {}
		
		# get start, end and strand of read
		start  = int(it[1])
		end    = int(it[2])
		strand = it[5]
		
		if extend == 0: extend = end - start ## no extending is needed, such as DNase-seq or FAIRE-seq
		
		weight = 1
		if (read_weight == True): weight = float(it[4]) ## the weight of a read ranging from 0 to 1
		
		# create tag by extending from 5'end of read
		tag_start = start
		tag_end   = start + extend
		if strand == "-":
			tag_start = end - extend
			tag_end   = end
		if tag_start < 0: tag_start = 0
		
		# Put back to dictionary
		for i in xrange(tag_start, tag_end):
			if not outdict[chr].has_key(i): outdict[chr][i] = weight
			else: outdict[chr][i] += weight

	return (outdict, totalcount)




def nuDict(bed):
	#print >> sys.stderr, "-Build nucleosome bed dictionary..."
	nucledict = {}
	if os.path.isfile(bed):
		bedDict (nucledict,path = bed, header = 0)
	else:
		bed = bed + "/"
		bedfiles = listDir(bed,"bed")
		for bedfile in bedfiles:	
			print bedfile
			bedDict (nucledict,path = bed + bedfile, header = 0)
	return nucledict

## Count the occupancy of a genomic factor such as H2AX around promoter
def promoter_occupancy (tss_path, bed_gf, winsize = 100, step = 20):
	outfile = open( os.path.basename(tss_path) + "." + os.path.basename(bed_gf) + ".count", 'w')
	
	print >>sys.stderr, "Loading GF dictionary..."
	dict_gf = twostrandDict(bed_gf)

	print >>sys.stderr, "Loading promoter dictionary..."
	dict_promoter = {}
	bedDict (dict_promoter, tss_path, header = 0)
	
	print >>sys.stderr, "Processing each chromosome..."
	
	for chr in chromList_new():
		if dict_gf.has_key(chr) and dict_promoter.has_key(chr):
			print >>sys.stderr, "Processing each promoter in", chr
			tss_sites = dict_promoter[chr]

			for site in tss_sites:
				
				items = site[0:-1].split("\t")				
				left  = int(items[1]) ## left pos of promoter
				right = int(items[2]) ## right pos of promoter
				gene = items[3]
				expr = items[4]
				
				outline = gene + "\t" + expr
				for i in range ( left, (right + step), step):
					start = i #start pos of window
					end   = i + winsize # end pos of window
					count = countReads(start, end, dict_gf[chr])
					outline = outline + "\t" + str(count)
				outline += "\n"
				outfile.write(outline)
	
	outfile.close()
	

def savedict (path, dict):
#This saves dictionary in a file
	file = open (path, 'w')
	pickle.dump(dict, file)
	file.close()

def nucle(bed = '',tss='',outpath='',windows=200,step=5,hd=1,scol=5,norm=1,argv = None,nucledict = {},fdict={},rdict={},writemode= 'w'):
#This function is used to calculate nucleosome occupancy around a genomic feature such as TSS, hotspot, donor, accepter.
#- Nucleosome is in BED file or a directory of BED file (one for a chromosome)
#e.g. 
#chr1    63475684        63475764        PSI179204_0002:3:1:1114:8170#0  0       +
#chr1    109778285       109778365       PSI179204_0002:3:1:1114:7473#0  4       -

#- Genomic feature is a file of chromsome coordinate of TSS, donor,.... This fuction considers the coordinate as center point then
#caculate nucleosome occupancy for left and right sides of the center, each <windows> number of windows, the width of window = step.
#e.g.
#chrom   txStart
#chr1    938709
#chr1    6596342

# fpath is path of saved dictionary of forward nucleosome
# rpath is path of saved dictionary of reverse nucleosome
	##~~~Handle input data
	if not argv is None:
		message = """
	This program gets the nucleosome score around genomic feature such as TSS, hotspot, exon
	Usage:
		<program>.py <Options> 
	Options:
		[-bed: nucle BED file path/string: ' ']
		[-pck: path of saved nucleosome dictionary: ' ', e.g. H2B_unique.pck]
		<-gf: genomic feature path/string> 
		< -o: output path/string>
		[ -w: num of windows: 200] 
		[ -s: step: 5] 
		[ -sc: strand column of nucle BED file: 5]
		[ -h: gf header: 1]
		[ -norm: normalization factor for each window: 1]
		"""

		para_dict = {'-bed':' ','-pck':' ','-gf':'','-o':'','-w':200, '-s':5,'-h':1,'-sc':5,'-norm':1}
		para_type = {'-bed':'str','-pck':'str','-gf':'str','-o':'str','-w':'int','-s':'int','-h':'int','-sc':'int','-norm':'int'}
		if para_process (argv[1:len(argv)], para_dict, para_type, message)==-1: sys.exit()

		bed = para_dict['-bed']
		tss = para_dict['-gf']
		outpath = para_dict['-o'] + "/"
		windows = para_dict['-w']
		step = para_dict['-s']
		hd = para_dict['-h']
		scol = para_dict['-sc']
		norm = para_dict['-norm']
		pck =  para_dict['-pck'] 
		fpath = pck + ".fw"
		rpath = pck + ".rv"
		##~~~~~~


	##MAIN~~~~
	if not os.path.exists(outpath): os.makedirs(outpath)
	fout = open (outpath + "forward.nuc", writemode) #out profile for reads mapped to forward strand
	rout = open (outpath + "reverse.nuc", writemode) #out profile for reads mapped to reverse strand
	
	#Build nucleosome dictionary
	if os.path.isfile (fpath) and os.path.isfile (rpath):
		print >> sys.stderr, "-Loading dictionaries from files..."
		fdict = pickle.load( open(fpath) )
		rdict = pickle.load( open(rpath) )

	elif nucledict == {} and bed != ' ':
		nucledict = nuDict(bed)
		[fdict,rdict] = strandDict(nucledict, col = 1, scol = scol)
	else:
		print >> sys.stderr, "Error: please check nucleosome path!"
		sys.exit()
	#Build feature dictionary
	print >> sys.stderr, "-Building genomic feature dictionary...."
	tssdict = {}
	bedDict (tssdict, path = tss, header = hd)
	print >> sys.stderr, "-Processing each chromosome...."
	for chr in tssdict.keys():
		if not fdict.has_key(chr): continue
		print >> sys.stderr, "  " + chr	
		for tssline in tssdict[chr]:
			tss_pos = int(tssline.split()[1])
			for win in range(-windows,windows):
				left = win*step + tss_pos
				right = (win+1)*step + tss_pos
				
				ftags = 0.0
				rtags = 0.0
				for i in range (left,right):
					if fdict[chr].has_key(i):
						ftags += fdict[chr][i]
					if rdict[chr].has_key(i):
						rtags += rdict[chr][i]
				ftags = ftags/norm #normalization by dividing window length = step
				rtags = rtags/norm
				fout.write (str(ftags) + '\t')
				rout.write (str(rtags) + '\t')
			fout.write ('\n')
			rout.write ('\n')

	fout.close()
	rout.close()


##----------------------------------------------
def windowsBed (chr, left, right, size, step, lastwin = False):
## lastwin is to decide whether to include the last window which has smaller size than the others
## step is sliding step, if it is equal to size --> moving window
	bed = []	
	wins = range(left, right, step)
	if not lastwin:
		wins = wins[:-1]
	k=0
	for i in wins:
		k += 1
		start = str(i)
		end   = str(i + size - 1)
		
		if lastwin and k==len(wins) : end = right + 1 ## added on 22-Jul-2011
		
		line  = chr + "\t" + start + "\t" + str(end) + "\n"
		bed.append(line)
	
	return bed


##----------------------------------------------
# Calculates proportion of bases covered by BED file
def get_Ndict(ambi_path): ##get number of ambiguous bases in each chrom
	ndict = {}
	lines = open(ambi_path).readlines()
	for l in lines:
		it = l.strip().split("\t")
		ndict[it[0]] = int(it[-1])
	return ndict
	
def coverage(lines, chr, chrlen, nlen, shift, sep):	
	dict = {}
	for line in lines:
		items = line.strip().split(sep)		
		
		if items[5] == "+":
			start = int(items[1])
			end   = int(items[1]) + shift
		elif items[5] == "-":
			start = int(items[2]) - shift
			end   = int(items[2])
		else:
			print >> sys.stderr, "Something wrong with strand!"
			sys.exit()
		
		if start < 0: start = 0
		if end > chrlen: end = chrlen
		
		for i in range( start, end ):
			if not dict.has_key(i): dict[i] = 1
			else: dict[i] += 1
		
	cover   = 100 * ( len(dict)/(float(chrlen)) )
	cover_N = 100 * ( len(dict)/(float(chrlen - nlen)) )
	print >> sys.stderr, chr + "\t" + str( round(cover,2) ) + "\t" + str( round(cover_N,2) )
	
def genome_coverage (bedfile, sep="\t", skip= 0, shift = 160, hg="hg18"):	
	if hg=="hg18":
		chrdict = get_chrominfo(path="/data2/thong/human_index/chromInfo_hg18.txt")
		ndict   = get_Ndict("/data2/thong/human_index/hg18.fa.ct")				
	else:
		chrdict = get_chrominfo(path="/data2/thong/human_index/chromInfo_hg19.txt")
		ndict   = get_Ndict("/data2/thong/human_index/hg19.fa.ct")		
	
	bdict = {}
	bedDict (bdict, bedfile, col=0, change_name = False, header = 0, sep = sep)
	if not bdict:
		print >> sys.stderr, "The file is empty!"
		sys.exit()
	for chr in chromList_new():
		if not bdict.has_key(chr): continue
		coverage(bdict[chr], chr, chrdict[chr], ndict[chr], shift, sep)

## Get stdin arguments from sys.argv[]
def arg (n, type = "str", exit = False, mes = "", default = ""):
	try:
		if type.lower().find("num") != -1:
			return float(sys.argv[n])
		return sys.argv[n]
	except:
		if exit:
			print >> sys.stderr, "Arg", n, "is wrong!", mes
			sys.exit()
		else:
			return default

## Re-write the function countWindow
## Count #reads belonging to each equal-size window cross genome (non-overlapping windows)
## Please use genome_density.sh for sliding windows.
## shift = 1/2 nucleosome size. skip: number of top lines to skip
def count_moving_windows(chromfile, bedfile, winsize, shift = 75, skip = 0, sep = '\t', lastwin = False):
	enrich = {}

	## Get chromosome sizes based on UCSC chromsome size information	
	chrdict = get_chrominfo(path = chromfile)
		
	chroms = chromList_new()

	## Get bed dict
	print >> sys.stderr, "Processing each tag from BED...."
	bed = getStream(bedfile, header = skip)

	#process each read (tag) in bed file
	for tag in bed:
		temp = tag[0:-1].split(sep) #remove \n and split
		chr = temp[0] #chromsome name in the first column of BED file
		
		if chroms.count(chr) == 0: continue #not in the chrom list 1, 2,..., 22, X, Y, M
		
		#initiation chr list (by a zeros vector) that stores the number of reads in each window of chr
		#the number of reads in a window is stored as an item of chr list
		if not enrich.has_key(chr):			
			numwin = chrdict[chr]/winsize + 1  #the last window may not have the size of 'winsize'
			if not lastwin: numwin = numwin - 1
			enrich[chr] = [ 0 for i in range(numwin) ]

		#find the center of nucleosome by moving 75bp forward or backward
		if temp[5] == '+': #forward strand
			center = int(temp[1]) + shift
		elif temp[5] == '-': #reverse strand
			center = int(temp[2]) - shift
		else:
			print >> sys.stderr, "Something wrong with strand!"
			sys.exit()
			
		win = center / winsize #get the index (window) of list that the tag belongs to
				
		if win < len(enrich[chr]):			
			enrich[chr][win] += 1

	return enrich
